Viscosity and viscosity anomalies of model silicates and magmas: a numerical investigation 
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We present results for transport properties (diffusion and viscosity) using computer simulations. Focus is 
made on a densified binary sodium disilicate 2Si02-Na20 (NS2) liquid and on multicomponent magmatic liq- 
uids (MORB, basalt). In the NS2 liquid, results show that a certain number of anomalies appear when the 
system is densified: the usual diffusivity maxima/minima is found for the network-forming ions (Si,0) whereas 
the sodium atom displays three distinct regimes for diffusion. Some of these features can be correlated with 
the obtained viscosity anomaly under pressure, the latter being be fairly well reproduced from the simulated 
diffusion constant. In model magmas (MORB liquid), we find a plateau followed by a continuous increase of 
the viscosity with pressure. Finally, having computed both diffusion and viscosity independently, we can dis- 
cuss the validity of the Eyring equation for viscosity which relates diffusion and viscosity. It is shown that it 
can be considered as valid in melts with a high viscosity. On the overall, these results highlight the difficulty of 
establishing a firm relationship between dynamics, structure and thermodynamics in complex liquids. 
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I. INTRODUCTION 

Viscosity is one of the key properties influencing the overall 
behavior of magmatic liquids. It is a fundamental property in 
Earth Sciences and has therefore led after nearly five decades 
of intensive research to a huge body of experimental and the- 
oretical studies and data base (Giordano et al., 2008; My sen 
and Richet, 2005). Systematic investigations with compo- 
sition and temperature in simple or complex silicates have 
been performed, and some generic trends have been iden- 
tified that have become quite popular. For instance, it has 
been found that numerous systems were displaying an Ar- 
rhenius behavior (Micoulaut, 2010) with temperature of the 
form r| = r|oexp[EA/A:BT] where Ea represents the activation 
energy for viscous flow. When properely rescaled with a ref- 
erence temperature, T g , at which the viscosity is 10 12 Poise, 
two categories of liquids have been identified (Angell, 1995) 
: those which remain Arrhenius-like in the super-cooled liq- 
uid over the entire range of temperatures implying that Ea 
does not depend on T, and those which display a curvature 
in the Arrhenius representation (a semi-log plot as a func- 
tion of inverse temperature) indicating that one has a tem- 
perature dependence in the activation energy for viscous flow 
(Wang et al., 2006). Usually, such liquids are fitted quite 
successfully with a Tamman-Vogel-Fulcher (TVF) law (Path- 
manathan and Johari, 1990; Giordano and Dingwell, 2003a,b) 
: r| = r|oexp(A/(T — To), A and To being constants. 

There has been a continuing effort to measure viscosity 
as function of pressure (Scarfe, 1973; Kushiro, 1976, 1978; 
1986; Kushiro et al., 1976; Scarfe et al., 1979; Brearley et 
al., 1986; Dunn and Scarfe, 1986; Mori et al., 2000; Suzuki 
et al., 2002, 2005, 2011; Reid et al., 2001; 2003; Tinker et 
al., 2004; Liebske et al., 2005; Ardia et al., 2008; Del Gaudio 



"Corresponding author: mmi@lptmc.jussieu.fr 



and Behrens, 2009). At constant temperature, the evolution of 
viscosity with pressure is complex and depends on the com- 
position of the silicate melt, in particular on its degree of de- 
polymerization. With highly polymerized melts (e.g. albite 
and dacite) the viscosity decreases continuously with pres- 
sure, though its behavior above the maximal pressure of in- 
vestigation (~7 GPa for albite and dacite) cannot be clearly 
anticipated. However, it is observed that the decrease of the 
viscosity with pressure becomes more gradual when the tem- 
perature is raised and can lead to a plateau value at high pres- 
sure (e.g. jadeite). When the melt is less polymerized (e.g. 
albite-diopside system), after an initial decrease, and accord- 
ing to the temperature, the viscosity exhibits a plateau value 
or a weakly pronounced minimum with the pressure. With 
depolymerized melts (e.g. diopside and peridotite) a different 
situation occurs. The viscosity first increases with rising pres- 
sure and reaches a maximum value prior to decrease at higher 
pressure, but higher the degree of depolymerization of the 
melt more pronounced is the viscosity maximum. Moreover, 
in a more general way, the pressure evolution of the viscosity 
has been found to be (anti-) correlated with that of the oxygen 
and silicon diffusion coefficients (Shimizu and Kushiro, 1991; 
Poe et al., 1997; Reid et al., 2001; Tinker et al., 2003), a find- 
ing which has been interpreted as the signature of a pressure- 
induced structural rearrangement. 

The understanding of changes in the viscous flow with ap- 
plied pressure can be suitably investigated by molecular dy- 
namics (MD) simulations. There are several ways to compute 
viscosity. Readers should refer to Allen and Tildesley (1987) 
for a detailed discussion and presentation. At equilibrium, the 
computation of viscosity from MD simulations can be per- 
formed by using the Green-Kubo (GK) formalism which is 
based on the calculation of the stress auto-correlation func- 
tion (Boon and Yip, 1980) given by: 

F(t) = (Pop(0^ap(0)> (1) 
where P a p(t ) is the ap component (a, p)=(x,y,z) of the molec- 
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ular stress tensor defined by: 

Pa^t^f + tt^tj «^P' ( 2 ) 

Ml Mly>! 

where F" is the component a of the force between the ions i 

and j, rf, and vf being the p component of the distance be- 
tween two atoms i and j, and the velocity of atom i, respec- 
tively. Alternative forms for P a p exist (Allen and Tildesley, 
1987) but these are found to lead to spurious effects when peri- 
odic boundary conditions are considered (Lee, 2007) or when 
the Ewald sum is not performed in a neat way (Alejandre et 
al., 1995; Allen et al., 1994). In this case, F(t) is found to not 
decay to zero in the long-time limit, as expected (Lee, 2007; 
Green, 1957; Kubo, 1957). One other possibility is to com- 
pute the viscosity within the framework of Non-Equilibrium 
Molecular Dynamics (Lees and Edwards, 1972). Here, the 
ratio of shear stress to a strain rate is computed and extrapo- 
lated to the limit of zero driving force (Borodin et al., 2009; 
Fernandez et al. 2004; Cherne and Deymier, 1998). 

Coming back to the GK formalism commonly used with 
MD simulations, the viscosity is calculated from the time in- 
tegral of the stress auto-correlation function F(t) (see eqs. (1) 
and (2)), 

The viscosity of a few silicate melts has been evaluated by 
MD simulation in following this route. Ogawa et al. (1990) 
were the first to evaluate the viscosity of a silicate melt (a 
sodium disilicate) by MD simulation but their results are un- 
certain due to a poor statistics. For silica, Horbach and Kob 
(1999) determined the temperature behavior and Barrat et al. 
(1997) studied the pressure dependence at high temperatures 
but without pointing out explicitly that their results showed a 
viscosity minimum at about 20 GPa. More recently Lacks et 
al. (2007) in investigating by MD the properties of silicate liq- 
uids along the MgO-Si02 join found a viscosity minimum at 
about 20 GPa in pure silica, a viscosity minimum which disap- 
pears progressively with increasing MgO content. Since then, 
these results have been confirmed by the MD studies of Ad- 
jaoud et al. (2008, 201 1) on Mg2Si04 melts, those of Spera et 
al. (201 1) on molten enstatite and also by first-principles MD 
calculations of Karki and Stixrude (2010a,b) on liquid silica 
and on liquid enstatite. We finally mention the recent studies 
of Karki et al. (201 1) and Gosh and Karki (201 1) on anorthite 
liquid and Mg2Si04 liquid, respectively. 

II. MODELING SILICATES 

A. Simulation details 

To complement the above MD studies and to bring some 
new information, we have investigated in the present study 
two silicate melts exhibiting comparable NBO/T ratio (~ 
0.8 and 1) but with very different chemical compositions; 



a sodium disilicate (NS2) and a mid-ocean ridge basalt 
(MORB). The molecular dynamics (MD) simulations were 
performed with the DL_POLY 2.0 code (Smith and Forrester, 
1996) in NPT Ensemble. The equations of motions for atoms 
were solved with a time step of 1-2 fs (10~ 15 s) by the leap- 
frog algorithm. 

Liquid Na20-2Si02 has been simulated by placing 666 sil- 
icon, 666 sodium, and 1668 oxygen atoms in a cubic box 
with periodic boundary conditions. As the atoms are bearing 
charges (ions), the long range coulombic forces are evaluated 
by a Ewald sum. The atoms interact via a two-body poten- 
tial (Born-Majer type) parametrized by Teter (2003, see also 
Cormack et al., 2003), which writes 

Vij(r) =A lj .exp(-r/p ii )-C ij /r 6 +z,z j /r, (4) 

where r is the distance between atoms i and j, zt is the effec- 
tive charge associated with the ion i, and where Ajj, p,j and 
Cij are parameters describing repulsive and dispersive forces 
between the ions i and j. In the present report, we do not focus 
on structure and thermodynamic properties of liquid NS2 sim- 
ulated with this potential and for a detailed analysis and com- 
parison with experimental data we refer the reader to Du and 
Cormack (2004, 2005). Notice that at zero pressure and 300K 
the density of the simulated NS2 glass is equal to 2.45 g/cm 3 , 
i.e. quite close to the experimental value (Vaills et al., 2001) 
2.37 g/cm 3 . 

The mid-ocean ridge basalt (MORB), a nine component 
system (0.13 wt% K 2 0, 2.94 wt% Na 2 0, 11.87 wt% CaO, 
7.77 wt% MgO, 8.39 wt% FeO, 1.15 wt% Fe 2 3 , 15.11 wt% 
A1 2 3 , 1.52 wt% Ti0 2 and 50.59 wt% Si0 2 ), has been sim- 
ulated with a two-body potential of the same form as eq. (4). 
The chemical composition of the MORB melt and the poten- 
tial parameters are those given by Guillot and Sator (2007a) 
in their simulation study of natural silicate melts (see Tables 
1 and 2 therein). The simulated sample of MORB was com- 
posed of 3,000 atoms (555 Si, 12 Ti, 195 Al, 9 Fe 3 +, 60 Fe 2+ , 
126 Mg, 141 Ca, 63 Na, 3 K, and 1,836 O) contained in a 
cubic box with periodic boundary conditions. The control pa- 
rameters (time step, Ensemble, etc.) of the simulation are the 
same as for liquid NS2. The density of MORB at liquidus 
temperature and its equation of state (i.e. density versus pres- 
sure) are well reproduced by the MD calculation (for a de- 
tailed description of the properties of the MORB melt, see 
Guillot and Sator 2007a,b). Notice that long simulation runs 
(~10 ns or 10 7 time steps) were performed to reach a good 
statistics when evaluating the viscosity by the Green-Kubo in- 
tegrand. 

Before presenting results, it is worthwhile to keep in mind 
some key data when one is about to calculate the viscosity of a 
supercooled melt. In using a high performance computer, the 
longest MD run that can be performed (in a reasonable time) 
for a system composed of ~ 3,000 atoms is about 1 /j s (or 10 9 
MD steps). To observe on this timescale the diffusion of atoms 
in the simulation box, the mean square displacement of each 
atom after 1 /j s has to be (on average) of the order of ~ 10 A 2 
(the mean distance between two nearest neighbors is roughly 
3 A in a melt, so the diffusive regime is reached when each 
atom has interchanged its position with a nearest neighbor). 
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Therefore, the smallest diffusion coefficient that can be eval- 
uated over a 1 /j s MD run is equal to ~ 2 x 10~ 14 m 2 /s (from 
Einstein relation D min = {r^ nin )/6t max where (rf nin )=lQ A 2 and 
imax = 10~ 6 s). In applying the Eyring theory to viscous flow 
, the viscosity r| can be deduced from the diffusion coefficient 
D of network forming ions through the equation, 



11 = 



kjT 
ID 



(5) 



where T is the temperature and A, a jump distance in the vis- 
cous flow. Experimentally it has been shown (Shimizu and 
Kushiro, 1984) that this equation works well with viscous liq- 
uids (e.g. silicates with high silica contents) provided that the 
D value is assigned to that of O or Si atoms and X ~2.8A, 
the oxygen-oxygen mean distance in silicate melts. Using 
equ. (5), a value of Do equal to 2xl0~ 14 m 2 /s leads to a 
corresponding viscosity af about 3000 Pa.s at 1473 K. We 
can use the Maxwell relation for which the relaxation time 
for viscous flow is related to the viscosity through the equa- 
tion, Xreiax = T|/Goo, where Goo ~ 10 10 Pa.s is the shear mod- 
ulus at infinite frequency (Dingwell and Webb, 1989). So to 
fully relax a viscous melt of viscosity r| ~ 10 4 Pa.s one needs 
to perform a MD run longer than X re i m = 1/js, a finding con- 
sistent with the evaluation based on the diffusion. Although 
qualitative, these estimates are useful because they point out 
that only the low viscosity regime (i.e. r| <C 10 4 Pa.s) can 
be investigated by MD simulation, the highly viscous regime 
(10 4 - 10 12 Pa.s) being unreachable with usual numerical re- 
sources. In the present study our MD runs did not exceed 
~10 ns, which means that only a viscosity value smaller than 
~10 Pa.s (100 Poise) could be investigated (notice that most 
of the MD studies published in the literature often do not ex- 
ceed 1 ns run intervals). 



B. Diffusion coefficients 



100 
10 

! 1 
0.1 
i 0.01 
0.001 
0.0001 

100 
10 

! 1 
0.1 
i 0.01 
0.001 
0.0001 




Na : Negodaev et al. 
Na : Gupta and King 
Na : Johnson et al. 
Si, O : Truhlarova et al. 
MD : Horbach et al. 
MD : Present work 



□ 







rprrrr 




10 12 14 16 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
O Si, O : Lesher et al. 
□ Na : Henderson et al. 
— MD : Present work 











1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

10 4 /T 



Figure 1: (Color online) A) Computed diffusion constants D^ fl , D& 
and Do in the NS2 liquid as a function of inverse temperature (blue 
curves and symbols), compared to experimental data for D^ a (Gupta 
and King, 1967; Negodaev et al. 1972; Johnson et al., 1951) and D S( , 
Do (Truhlarova et al., 1970), and to the simulated values of Dm,, 
Dj; and Do using an alternative potential (red curves and symbols, 
Horbach et al., 2001). B) Computed diffusion constants Djv fl . D^, 
and Do in the MORB liquid as a function of inverse temperature 
(blue curves and symbols), compared to experimental data for D^ a 
(Henderson et al.) and Dj,-, Dq (Lesher et al. 1996). 



We first compute the mean-square displacement of a tagged 
atom of type a in the melt, given by 

(^(OH^BMO-rKo)! 2 ) , (6) 

and extract from the dependence of (r 2 (t)} the long time 
behavior where the dynamics becomes diffusive. Using the 
Einstein relation lim r ^oo(r 2 (f))/6f = D, one can indeed have 
access to the diffusion constants D, (;=Si, O, Na) from the 
mean square displacement. These quantities are plotted for 
NS2 and MORB in Fig. 1 as a function of inverse temper- 
ature. Because of the slowing down of the dynamics at low 
temperature, the computation of D,- is here restricted to T > 
1500 K with NS2 and T > 1 173 K with MORB. For NS2, the 
diffusion of sodium atoms is found to be remarkably close to 
the experimental data obtained by Gupta and King (1967) and 
Negodaev et al. (1972), with a very good agreement some- 
where around IQ 4 /T ~ 6. For network forming ions (Si,0) 
an excellent agreement is found with the experimental data 



of Truhlarova et al. (1970) who determined the (Si,0) dif- 
fusion from a Si02 dissolvation reaction in a NS2 liquid. In 
this respect, the Teter potential appears to be very accurate for 
reproducing the diffusion of elements in the NS2 melt. We re- 
mark that the diffusion coefficients display an Arrhenius-like 
dependence when 10 4 /T>4 with activation energies equal to 
1.18 eV, 1.17 eV and 0.43 eV for Si, O and Na respectively. 

To be complete, one should stress that alternative poten- 
tials (Kramer et al., 1991) lead to diffusion coefficients for 
(Na,Si,0) that are found (Horbach et al., 2001) to be at least 
one order of magnitude lower than the experimental values 
(for 10 4 /T>5), and our own results on the MORB liquid as 
well (see below), whereas the corresponding activation ener- 
gies also differ by a factor of ~2 (Fig. 1). These findings 
illustrate how difficult the numerical reproduction of diffu- 
sion coefficients in silicate melts can be. This well-known 
feature has been highlighted for the case of silica by Hemmati 
and Angell (1997) who showed that diffusion is highly model 
dependent with differences between models increasing up to 
at least three orders of magnitude when the liquid is cooled 
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Figure 2: (Color online) Stress auto-correlation function F(t) of a 
2000 K liquid NS2 as a function of time for p=2.5 g/cm 3 , and a 2273 
K MORB liquid at the same density. 



down. So, the very good agreement between experimental 
and simulated ionic diffusivities in the NS2 melt, allows one 
to be reasonably confident in the ability of the Teter potential 
to describe accurately other transport properties as the melt 
viscosity. As for the MORB melt, the temperature depen- 
dence of the O, Si and Na diffusion coefficients is shown in 
panel B of Fig. 1 (among all the network modifying cations in 
MORB, only the diffusion coefficient of Na is shown because 
it is the most mobile ion and because it is common with the 
NS2 melt). Na is found to diffuse more slowly in MORB than 
in NS2 whereas Si and O atoms are diffusing faster. However, 
the agreement between simulated and experimental values is 
poorer with MORB than with NS2 (see Fig. 1 for a compar- 
ison). Actually, Na, Si and O atoms diffuse too rapidly in 
the simulated MORB (by a factor of ~5 for Na and ~50 for 
Si and O ). So one may anticipate some deviation between 
calculated and experimental viscosities for the basaltic liquid, 
and a better agreement for liquid NS2. 



C. Viscosity 

In Figure 2, we represent the stress auto-correlation func- 
tion F(t) of NS2 and MORB at around zero pressure, com- 
puted following equs. (1) and (2). As mentioned above, 
the stress auto-correlation function is indeed found to decay 
to zero. Typical oscillations which are found at the subpi- 
cosecond timescale are associated with frequencies in the Ter- 
ahertz range and above [10-300 THz] corresponding to the 
main band found in the vibrational density of states. This 
band is attributed to stretching motions of Si-O and Na-O 
bonds (and other cation-oxygen bonds in MORB). Note that 
these oscillations are damped as the pressure is increased (not 
shown). Following eq. (3), the time integration of the stress 
auto-correlation function leads to the viscosity (Fig. 3). The 
viscosity of NS2 exhibits an Arrhenius-like behavior over a 
broad temperature range (2<10 4 /T<6). The calculated acti- 
vation energy for viscous flow is about ~ 1 .33 eV which is in a 



reasonable agreement with Bockris data (~ 1 .65 eV) measured 
(Bockris et al., 1954) in the temperature range (1723-1373 K 
or 5.8<10 4 /T<7.3). Furthermore, in the temperature range 
where calculated and experimental values can be compared 
with each other (i.e. 5.8<10 4 /T<6.7) calculated values are 
smaller by a factor of ~2 than the experimental data of Bock- 
ris et al. (1955). This deviation is surpring given the good 
agreement obtained for the diffusion coefficients (see Figure 
1), and given the fact that r| and D can be simply related (equ. 
5). 

In the case of MORB the temperature dependence of the 
calculated viscosity is found to be non Arrhenian and can 
be fitted accurately over a large temperature range by the 
TVF equation with r)oo=6.5x 10~ 3 Poise, A = 0.42 eV, and 
7b = 670 K. Experimentally, the viscosity of basaltic melts 
is known to be non Arrhenian, so the simulation reproduces 
this feature, but the simulated MORB is much less viscous 
than the real one (a factor of 50 at 1673 K, see Fig. 3), this 
is why A and To values are slightly different from those rec- 
ommended for basalts by the regression formula of Giordano 
and Dingwell (2003a,b), and Giordano et al. (2008). In fact 
the MORB melt simulated at 1273 K corresponds, as far as 
the viscosity is concerned, to the real melt at about 1673K, 
this temperature deviation diminishing at higher temperatures 
where the viscosity of any silicate melt (whatever is its com- 
position) tends to be very low (less than 0. 1 Poise above 3000 
K). More generally, the agreement between experimental and 
calculated viscosities for silicate melts is actually difficult to 
obtain because the viscosity seems to be even more sensitive 
to the details of the potential than the diffusion coefficients 
themselves (for a related discussion see Vuilleumier et al., 
2009). Moreover, direct comparisons between experimental 
and simulated values of the viscosity are scarce in the litera- 
ture (e.g. Lacks et al., 2007; Adjaoud et al., 2008; Karki and 
Stixrude, 2010b, Spera et al., 2011) and are usually done at 
very high temperature (>2500K) where the viscosity is very 
low and its experimental value badly known or obtained from 
uncontrolled high-temperature extrapolations. 

In the TVF equation, rioo represents the high temperature 
limit to silicate melt viscosity. By analyzing measured viscos- 
ity curves for a great number of silicate liquids, it has been 
shown that all silicate melts converge to a common value of 
viscosity at very high-T, about 10 2 -10 3 Poise (Russell et 
al., 2003; Giordano et al., 2008, Zheng et al., 2011). But it 
is difficult to test this conjecture because of a lack of viscos- 
ity data for silicate liquids at super-liquidus temperatures (i.e. 
for T > 2000 K). In contrast, the simulation can shed some 
light on this topic because the high-T range is easily acces- 
sible in a numerical experiment. However it is worthwhile 
to clarify what means the concept of a high-temperature limit 
for a silicate liquid. When a silicate melt is heated isobarically 
(e.g. under atmospheric pressure) it vaporizes incongruently 
in forming a gaseous phase which is a complex mixture of 
atoms and molecules (e.g. the major vapor species above liq- 
uid silica are SiO, Si02, O2 and O, and about thirty species 
compose the vapor phase above a basaltic lava, for a review 
see Schaefer and Fegley, 2004). Nevertheless this is not a 
real problem as long as the resulting vapor pressure remains 
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low, i.e. the temperature is not too high (e.g. the vapor pres- 
sure above basaltic lavas is 10~ 5 -10~ 4 bar at 1900 K). This 
situation changes drastically at higher temperature where the 
vapor pressure becomes high (e.g. ~ 10~ 2 bar at 2500K above 
for a tholeiitic liquid). Then the vaporized fraction becomes 
rapidly non-negligible (10% or more) and the composition of 
the residual melt deviates significantly from the original one. 
To avoid this problem, one only has to prevent vapor for- 
mation in applying on the liquid an isostatic pressure larger 
than the saturating vapor pressure. A convenient thermody- 
namic path consists to follow (or to be as close as possible) 
the liquid branch of the liquid-vapor coexistence curve up to 
the critical point. The critical point parameters for silicates 
have never been measured because of the very high tempera- 
tures involved (much beyond 5000 K, estimations for SiC>2 are 
given in Guissani and Guillot (1994) and Melosh (2007)). By 
imposing the pressure, we have been able to evaluate the vis- 
cosity of our simulated MORB along (or close by) the liquid 
branch of the liquid-vapor coexistence curve. However due 
to the density fluctuations generated by the MD calculation it 
becomes more and more difficult to preserve the cohesion of 
the liquid phase in rising the temperature because the liquid 
then tends to transform spontaneously into the corresponding 
vapor phase. So the highest temperature reached in keeping 
the cohesion of the liquid phase around 1 bar was 4673K, a 
temperature at which the density is equal to 1.47 g/cm 3 and 
the viscosity is evaluated to 10~ 28 Poise. This latter value 
matches rather well the estimated range for r|„o (10-2-10" 3 
Poise, see above). It is certainly possible to pursue the calcu- 
lation at higher temperature up to the critical point, but it is a 
tedious task to locate accurately the critical point of our model 
and we leave that to a future work. In using the same proce- 
dure with liquid NS2 we have reached 4500K, a temperature 
at which the liquid density is as low as 1 .66 g/cm 3 and its vis- 
cosity 10~ 25 Poise. We anticipate that the common value of 
T|oo for MORB and NS2 at the critical point could be about 
10~ 3 Poise or slightly less. Notice that in approaching the 
critical point, the structure of a silicate liquid is quite different 
of the structure from the melt near the liquidus. The near crit- 
ical liquid is indeed fully depolymerized and is composed of 
molecular clusters involving various species (oxides). So it is 
expected that all silicate liquids will exhibit similar flow be- 
havior in the critical region, that of a regular molecular liquid. 

We now follow both viscosity and diffusion with pressure 
in our two melts (NS2 and MORB). In Fig. 4 and 5, we plot 
the behavior of the diffusion constants with density and pres- 
sure. The most mobile atom, Na, is found to display the same 
behavior in the two melts: its diffusion coefficient is contin- 
uously decreasing with the pressure. In the case of MORB 
the other network modifying cations (Ca, K, Mg,...) show the 
same behavior with the pressure as Na does (Fig. 5). 

The behavior of network forming ions (Si,0) is quite dif- 
ferent as their diffusion coefficients go through a maximum 
value with the pressure. Furthermore this maximum is lo- 
cated at about 15 GPa in NS2 and about 50 GPa in MORB 
and it is much more pronounced in the NS2 melt than in the 
basaltic composition. Notice that a diffusivity maximum for 
network forming ions is also observed with silica (Shell et al. 
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Figure 3: (Color online) Simulated viscosity of the NS2 liquid at zero 
pressure (filled red squares), compared to experimental data (red cir- 
cles) from Bockris et al. (1955), together with simulated viscosity 
of the MORB liquid at zero pressure (filled black squares), com- 
pared to experimental data (broken black curve) from Villeneuve et 
al. (2008). The dotted line is a high temperature Arrhenius fit for the 
MORB data whereas the solid line is a TVF fit (see text for parame- 
ters). 
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Figure 4: (Color online) Silicon (red), Oxygen (green), and sodium 
(black) diffusion constants of a NS2 liquid at T=2000 K as a function 
of system density (bottom axis) and pressure (top axis). Right axis: 
Isothermal compressibility of the liquid computed from the equation 
of state at 2000 K. 



2002), germania (Dreyfus and Micoulaut, 2012), some sodo- 
aluminosilicates (Poe et al., 1997) and in water (Errington and 
Debenedetti, 2001). 

Is such an anomaly also visible in the evolution of viscosity 
with pressure ? Actually, for NS2 we do find that r| displays a 
minimum of viscosity but in a pressure range 4-10 GPa which 
does not coincide with the diffusivity maximum for Si and 
O (Pmtu=15 GPa). For MORB, the viscosity shows no mini- 
mum, it is nearly constant up to 5 GPa and increases progres- 
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Figure 5: (Color online) Silicon (red), Oxygen (green), and sodium 
(black) diffusion constants of a MORB liquid at T=2273 K as a func- 
tion of pressure. 
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Figure 6: (Color online) Simulated viscosity of a NS2 liquid at 
T=2000 K (filled circles) and a MORB liquid at T=2273 K (open 
circles) as a function of pressure. 

sively at a higher pressure. These simple findings suggest that 
the diffusivities of network forming ions are not the only ones 
responsible of the viscous flow, network modifying cations 
play also a role. 



III. DISCUSSION 

As briefly recalled in the introduction, the pressure de- 
pendence of the viscosity has been reported experimentally 
for a number of silicates (Scarfe, 1973; Kushiro, 1976, 
1978; 1986; Kushiro et al., 1976; Brearley et al., 1986; 
Dunn and Scarfe, 1986; Mori et al., 2000; Suzuki et al., 
2002, 2005, 2011; Reid et al., 2003; Tinker et al., 2004; 
Liebske et al., 2005; Ardia et al., 2008). The analysis 



of these data along the sequence of composition albite—^ 
jadeite— > rhyolite— > andesite^ basalts^ albite x -diopsidei_ x 
— > jadeite x -diopsidei_ v — > diopside — » peridotite, which cov- 
ers a large NBO/T range (from to 2.8), leads to a con- 
trasted picture. In the initial pressure range (0-2 GPa), the 
slope (dri/dP)^ is negative for composition with NBO/T < 
1.2 and is positive for higher NBO/T value. Although the 
NBO/T ratio is a good indicator of melt polymerization, its 
use for predicting the pressure evolution of viscosity requires 
some caution (e.g. Toplis and Dingwell, 2004). For instance, 
albite and jadeite melts are both considered as fully polymer- 
ized (NBO/T=0) but the pressure behavior of their viscosity 
is different (a viscosity curve presenting a weak minimum lo- 
cated about 2-3 GPa for jadeite as compared with a continuous 
decrease of the viscosity over several decades up to 7 GPa for 
albite, see Suzuki et al. 2002, 201 1). In keeping in mind this 
information, and according to the above data analysis, the vis- 
cosity of NS2 and MORB composition (NBO/T=l and 0.8) 
should exhibit a negative slope. For instance, Scarfe et al. 
(1987) notice that the viscosity of the NS2 melt is decreas- 
ing between and 2 GPa. In the same way, the viscosity of 
oceanic basalts (Kushiro et al., 1976; Kushiro, 1986; Ando, 
2003) maintained at fixed temperature decreases slightly with 
increasing pressure up to 3 GPa where it reaches a minimum 
value prior increase upon further compression. But as an in- 
crease of the temperature tends to decrease the viscosity of 
a melt maintained at a fixed pressure, the weak minimum of 
viscosity in basalts is expected to disappear in increasing the 
temperature. This is why a viscosity minimum is barely vis- 
ible with our MORB melt simulated at 2273K, especially as 
the latter one describes the viscous flow of a real MORB at a 
higher temperature, may be as high as 2600-2700 K (see Fig. 
3). In contrast the viscosity of the simulated NS2 melt at 2000 
K presents a clear minimum as a function of pressure (notice 
that the temperature dependence of the viscosity of the NS2 
melt is much better reproduced by the simulation than in the 
MORB case, see Fig. 3). 

Since the first report of an anomalous pressure dependence 
of the viscosity ((dri/dP)^ <0) by Kushiro (1976) and Kushiro 
et al. (1976), different authors have tried to explain the ori- 
gin of this phenomenon. Using the Adam-Gibbs theory and 
a pressure dependence for the degree of melt polymeriza- 
tion, Bottinga and Richet (1995) have shown that the decrease 
of the specific volume with alkali content was leading to an 
increased sensitivity of the viscosity with pressure change. 
Gupta (1987) has proposed that a system having a thermal 
expansion coefficient in a glass larger than that of a liquid 
could display such anomalies. However, this does not seem to 
be verified experimentally (Knoche et al., 1994). A relation- 
ship with structural changes under pressure also has been pro- 
posed, in particular a pressure-induced coordination change of 
Si and Al or Si-O and Al-O bond weakening due to pressure- 
induced distortion of the liquid structure (Waff, 1975; My sen 
et al., 1980). However, a study (Sharma et al., 1979) on dense 
liquid Ge02 where both coordination change (Micoulaut et al. 
2006) and a pressure induced viscosity anomaly occur could 
not establish the correlation. This does not rule out the pos- 
sibility of a structural origin for the anomaly. For instance, 
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Suzuki et al. (2002) do find a correlation between the gradual 
decrease of the viscosity of albite with pressure up to 7 GPa 
and the increasing concentration of 5 -coordinated Al. 

The variation of r| along an isotherm also suggests that the 
corresponding activation energy for viscous flow EA(P,T) is 
minimum. Here, for selected pressures in the NS2 system, we 
have computed r| with temperature and find that the activa- 
tion for viscous flow should minimize indeed in the pressure 
window, i.e. we obtain Ea=1.33, 1.09 and 1.34 eV for P=0, 
6, and 22 GPa, respectively. The link between minima 
and the presence of a stress-free state in corresponding low- 
temparature glassy networks has been established recently 
(Micoulaut, 2010), indicative of a flexible to rigid transition 
(Micoulaut and Phillips, 2003). Since it is known that at am- 
bient pressure the NS2 is flexible (Vaills et al., 2005), one may 
obtain a pressure-induced flexible to rigid transition with in- 
creasing P (Trachenko et al. 2004), the minimum in Ea (and 
r|) being one of the expected salient features. 

In summary, all these studies point out the subtle interplay 
between structure, dynamics and thermodynamics of the sili- 
cate melt that leads to anomalous behavior in transport prop- 
erties. 



A. Viscosity versus diffusion 

The connection between diffusion and viscosity can be re- 
alized in a simple way via the Eyring relationship given by eq. 
(5) involving X which is a typical jump distance for the diffus- 
ing atom. According to Eyring, the relationship (5) holds (i.e. 
X is constant and equal to the average interatomic distance) 
if the activated process for diffusion can be assumed identical 
with that of viscous flow. So for silicate melts the network 
forming ions (Si and O) are good candidates to fulfill this pre- 
requisite. In fact the Eyring relationship has been tested over a 
large range of melt composition, from simple binary oxides to 
natural magma compositions (Oishi et al., 1975; Yinnon and 
Cooper, 1980; Dunn, 1982; Shimizu and Kushiro, 1984; Dunn 
and Scarfe, 1986; Lesher et al., 1996; Reid et al., 2003; Tin- 
ker et al., 2004). In polymerized melts of high viscosity (e.g. 
jadeite), it is observed that the diffusivity of oxygen ions and 
the melt viscosity are related to each other through the Eyring 
relation in introducing for the jump distance, X, the average 
oxygen-oxygen distance in the melt (~ 2.8A). In depolymer- 
ized melts where the viscosity is generally lower (e.g. basalts 
and diopside), the Eyring relation is fulfilled with X larger than 
2-3 times the mean oxygen-oxygen distance. However even in 
polymerized melts, the Eyring relation fails to fully reproduce 
the correlated evolution of viscosity and oxygen self-diffusion 
coefficient. As a matter of fact it has been reported that oxy- 
gen self-diffusion coefficient in albite (Poe et al., 1997) and 
in dacite melts (Tinker and Lesher, 2001) reaches a maximum 
value at about 5 GPa whereas the viscosity of the correspond- 
ing liquids is continuously decreasing with the pressure and 
shows no viscosity minimum at 5 GPa (Suzuki et al., 2002; 
Tinker et al., 2004). In this context, we have evaluated the 
length, X = ksT /r\Do.si along two thermodynamic paths : 
the isobar P=0 (with varying temperature) and the isotherm 
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Figure 7: (Color online) Typical lengthscale X associated with the 
Eyring relation (5) as a function of the melt viscosity of a NS2 liq- 
uid, computed either from the pressure variation (at 2000 K, black 
symbols) or from the temperature variation (at P=0, blue symbols). 



T=2000K for NS2 (with varying pressure). For convenience 
this length is represented in Fig. 7 as function of the calcu- 
lated melt viscosity. Along each path (isobar or isotherm) the 
arrow indicates the direction of the increasing thermodynamic 
parameter (T or P). Notice that in the representation A.(r|) (Fig. 
7), both curves exhibit a maximum value (22A) for both the 
isotherm and the isobar which expresses a non-trivial behavior 
of r| with diffusivity. Clearly, for the two melt compositions 
and whatever the thermodynamic path, X decreases when the 
melt viscosity increases. A simple extrapolation suggests that 
the Eyring relation should be verified only when the viscosity 
is sufficiently high (for r| > 100 Poise), i.e. when A, becomes 
of the order of ~ 5A.. This conclusion is in agreement with 
the viscosity data discussed above. On the other hand, when 
the viscosity is becoming very low and the diffusion coeffi- 
cients of network forming ions high (> 10~ 9 m 2 /s), the length 
X matches better the values predicted by the Stokes-Einstein 
equation, X = 2nd, where d is the diameter of the diffusing 
molecule (~ 3 A for oxygen). As the Stokes-Einstein relation 
works well with molecular liquids (Li and Chang, 1955) it is 
not too surprising that it also leads to reasonable values for 
silicate liquids at very high temperature (e.g. in approaching 
the critical point) where the melt is essentially depolymerized 
and shares some structural similarities with regular liquids. In 
summary, the Eyring relation is not adapted to the description 
of the transport properties of silicate melts at superliquidus 
temperatures. Alternative models have been proposed in the 
literature (see Mungall, 2002) but a quantitative theory is still 
lacking. 



B. Structural correlations 

As mentioned before, there have been many attempts to cor- 
relate structural features with the viscosity anomalies in sili- 
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cate liquids (Waff, 1975; Woodcock and Angell, 1976; Mysen 
et al., 1980, Suzuki et al., 2002). From the liquid structures 
obtained for the two melt compositions (NS2 and MORB), we 
have computed the fraction of silicon atoms that are 4-, 5- and 
6-fold coordinated. Results are displayed in Fig. 8a. Here 
one sees that the melt structure evolves from roughly 100% 
4-coordinated Si™ to a growing proportion of Si v and Si vl , 
Si v increasing substantially for P larger than 10 GPa and go- 
ing through a maximum found at about 25 GPa in NS2 (and 
above this pressure for MORB) whereas Si 1 ' 7 is becoming sig- 
nificant only above 15GPa. 

From the probability distribution {p,} of Si' species 
(i=IV,V,VI), one can compute a configurational entropy ac- 
cording to: 



100 



-J^Pilapi, 



(7) 



which maximizes when the population of S v is maximum 
(~30 GPa). This configurational entropy is represented in Fig. 
8b as exp[l/S c ]. In fact, if one assumes that the Adam-Gibbs 
relation for viscosity of the form: r| = rjo exp[A/TS c ] is valid, 
one may expect from a simple inspection of the relationship 
that a minimum in r| is correlated with a maximum in S c , as 
also suggested by Bottinga and Richet (1995). However, for 
NS2 a quick look at figure 8b clearly shows that r|(P) and 
exp[l/S 6 ] are not correlated as the location of the maximum in 
S c does not coincide with the pressure range where a viscosity 
minimum is observed. In the case of MORB, the gradual in- 
crease of the viscosity with pressure could be correlated with 
that of S c . So it is hard to firmly conclude that a structural 
reorganization through coordination changes of Si is directly 
responsible of the pressure evolution of the viscosity. To be 
complete, notice that Goel et al. (2011) have recently shown 
by MD simulations that the diffusion coefficient of Si atoms in 
silicate melts along the MgO-Si02 join, correlates well with 
the excess entropy when the latter one is evaluated from the 
knowledge of the pair distribution functions. But it is unclear 
if this correlation also holds with the viscosity as the relation- 
ship linking diffusivity and viscosity is far to be obvious (see 
the above discussion). 

We have therefore calculated an excess entropy S2 in us- 
ing the expression of the two-body term (Baranyai and Evans, 
1989), which writes: 



S2 



k B p 



L 



<r/r 



(8) 



p being the density, and x, the concentration of species i. Re- 
sults of this calculation are shown in Fig. 8b (right axis). S2 
is found to decrease continuously with pressure and a corre- 
sponding Adam-Gibbs r| does not show a viscosity minimum. 
This behavior is at variance with that obtained by Goel et al. 
(201 1) in their MD simulation study of MgO-2Si02 (a system 
related to NS2) where S2 is exhibiting an entropy maximum 
around 10 GPa, very well correlated with a diffusivity max- 
imum, as also found by Jabes et al. (2010) for liquid Si02, 
Ge02, BeFe2 and H2O. So our search for an eventual relation- 
ship between viscosity and configurational entropy at the sin- 
gle temperature 2000 K is inconclusive and certainly deserves 
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Figure 8: (Color online) a) Fraction of 4-, 5- and 6-fold silicon atoms 
as a function of pressure in the NS2(black curves) and MORB liquid 
(red curves), b) Simulated viscosity of the NS2 liquid at T=2000 K 
(filled black circles, same as Fig. 3), compared to an Adam-Gibbs 
viscosity computed from a structure-related configurational entropy 
S c . Right axis: Entropy S2 (broken line) computed from equ. (8). c) 
Mean angle around the bridging oxygen (BO) atom and correspond- 
ing second moment (standard deviation) of the bond angle distribu- 
tion (right axis). 



a further scrutiny. An alternative way to understand the trends 
of the viscosity with pressure may arise from the evolution 
of the T-O-T intertetrahedral bond angle distribution (BAD) 
as it is well known that bond angles are much more pressure 
sensitive than bonds themselves (Micoulaut, 2004). Figure 
8c shows the evolution of the average value of the BAD be- 
tween a bridging oxygen (BO) and its two SKD4/2 tetrahedra, 
together with the second moment (standard deviation) of the 
BAD. The latter one quantifies the angular excursion around a 
mean value and it is expected to be pressure sensitive, as high- 
lighted in the case of densified GeSe2 (Bauchy et al., 2011). 
When the average value of the BAD decreases steadily in NS2 
melt from 150° to 140" between P=0 and 20 GPa, it varies be- 



9 



tween 140° and 128° in MORB over the same pressure range. 
In the same way the angular excursion increases (with pres- 
sure) equally well in the NS2 and in MORB melt as shown 
by the standard deviation Gq evolution with pressure. So, here 
again, it is not clear that viscosity and diffusion anomalies 
(e.g. in NS2) induced by the pressure can be univoquely asso- 
ciated with a particular structural rearrangement even if there 
is no doubt that structural changes are taking place in the sili- 
cate melt under investigation. 



IV. SUMMARY AND CONCLUSION 

In summary, we have used Molecular Dynamics simula- 
tions to study the dynamics of two systems of importance in 
materials science and geology : the sodium disilicate (NS2) 
and a basaltic (MORB) liquid. Results show that the diffu- 
sion constant computed from a two-body Born-Majer type 
potential at ambient conditions exhibits a remarkable agree- 
ment with experimental data for all species (Na, Si, O) in the 
NS2 liquid. This contrasts with the level of agreement ob- 
tained for the MORB system which shows a deviation by at 
least one order of magnitude when compared to experimen- 
tal findings. The same situation is nearly reproduced in the 
comparison of calculated viscosity using the Green-Kubo for- 
malism (based on the stress auto-correlation function) with 
experimental data from Bockris et al; (1955) for NS2, and 
Villeneuve et al. (2008) for MORB. 

We have then studied the behavior of both the diffusion con- 
stant and the viscosity with changing pressure in the two liq- 
uids with temperature close to 2000 K. In the NS2 system, a 
diffusivity maximum for the network forming species (Si,0) 
is obtained at around 15 GPa, a result which parallels similar 
findings for other densified tetrahedral liquids such as water, 
germania, or silica. The sodium cation is found to display 
three distinctive regimes for diffusion, the one corresponding 
to lower densities (negative pressure) being identified with a 
stretched melt (Bauchy and Micoulaut, 201 1). These features 
are not found in the more depolymerized MORB liquid which 
shows for Si and O diffusion a steady behavior between and 
5 GPa, followed by a continuous decrease applied pressure. 
Interestingly, a pressure window is found in the NS2 liquid 
with a minimum at around 5 GPa that does not seem to be cor- 
related with the maximum in the obtained diffusivity. Again, 
this contrasts with the steadily increase of the MORB liquid 
viscosity with pressure. 

Having in hand both the computed diffusion and the vis- 
cosity, we have investigated the validity of the Eyring equa- 
tion which relates both quantities and involves also a jump 
distance X usually of the order of the average oxygen-oxygen 
bond distance (2.8 A), as exemplified by numerous experi- 
mental studies. This investigation has been performed along 
two thermodynamic paths: at fixed pressure (P=0) and chang- 
ing temperature, and at fixed temperature (T~ 2000 K) and 
changing pressure. We have shown that A. calculated from the 
simulated diffusion and viscosity was indeed of the order of a 
few Angstroems, but only in the high viscosity regime (here 
100 Poise). It underscores the fact that the Eyring relationship 



should hold only for deeply supercooled liquids as it fails for 
liquidus temperatures. 

Finally, we have attempted to relate the obtained anomalies 
in viscosity with structural changes of the melt. With growing 
pressure, the population of higher coordinated silicon (Si , 
Si VI ) increases in a similar fashion in the NS2 and MORB 
liquid which rules out the possibility of a direct relationship 
between e.g. coordination increase and viscosity minimum. 
Moreover, in the NS2 liquid configurational entropies have 
been calculated and their behavior with pressure does not fol- 
low at all the trend obtained for the viscosity. More subtle 
structural and thermodynamic changes are clearly at play and 
one does not recover for NS2 the reported correlation between 
configurational entropy and diffusivity maxima found by Goel 
et al. (201 1) in the MgO-SiO? system. 

Our attempts to find out a quantitative link between 
pressure-induced structural rearrangement and viscosity 
anomalies is rather deceptive. Although some correlations be- 
tween diffusivity, viscosity and excess entropy certainly do 
exist and are sometimes emphasized in the literature, their ef- 
fectiveness depends strongly on the melt composition, a fact 
which makes risky a clear conclusion. On the other hand, a 
recent theoretical analysis by Schmelzer et al. (2005) seems 
to be promising. In this study the pressure dependence of the 
viscosity is analyzed from a thermodynamic al point of view. 
The working equation expresses the derivative of the viscosity 
with respect to the pressure under the following form, 

\dP) T a T \dT ) p\d^) T \dP) T ' W 

where i; is an order parameter (e.g. the degree of polymer- 
ization given by BO/(BO+NBO) where BO and NBO are the 
numbers of bridging and nonbridging oxygen, respectively), 
Kj is the compressibility and 0C7- the thermal expansion coef- 
ficient. If the order parameter is invariant with pressure (so 
d^/dP=0) then dr|/dP is governed by the first term in eq. (9) 
which expresses a free volume effect (r| only varies with the 
total volume of the melt). But as Kj is negative in a liquid 
at equilibrium, and dr\dT is positive (the viscosity decreases 
when the temperature increases), the variation of r| with P will 
be positive or negative according to the sign of the thermal ex- 
pansion coefficient, a^. If the latter one is negative then the 
viscosity will tend to decrease in increasing pressure and if a 7- 
becomes positive after a further compression then the viscos- 
ity will go to a minimum value before to increase. As shown 
by Schmelzer et al. (2005) this behavior is encountered with 
the viscosity of liquid water. However the first term in eq. (9) 
is not sufficient to reproduce all the magnitude of the viscosity 
minimum and the contribution of structural changes induced 
by the pressure has to be accounted for (via the second term in 
eq. (9)). In the case of silicate melts it is well documented that 
the viscosity increases when the degree of polymerization of 
the melt increases (so dr|/d^ >0) wheras for a given melt its 
degree of polymerization decreases when pressure increases 
(so d^/dP<0). Therefore the second term of eq. (9) is gen- 
erally negative. For NS2 and MORB we have checked that 
O.J is positive at liquid temperature and over a large pressure 
range which implies that the first term in eq. (9) is positive for 
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these two compositions. As our MD results for NS2 lead to a 
decrease of the viscosity between and ~4GPa, this suggests 
that the second term in eq. (9) overbalances the first term. 
For MORB the viscosity is found virtually constant in the 0- 
4 GPa range so the two terms nearly compensate each other 
in eq. (9). At high pressures (e.g. above 10 or 20 GPa) the 
two melts are more and more depolymerized and d^/dP be- 
comes small because the degree of depolymerization does not 
evolve very much with the pressure due to melt densification. 
So the free volume effect (first term) is the dominant contribu- 
tion to eq. (9) at high compression rates and the viscosity of 
the two melts increases steadily in the high pressure range. In 
summary it appears that eq. (9) is very useful to decipher the 
complex behavior exhibited by the viscosity of silicate melts 
as function of pressure. To be more quantitative and predictive 
it remains to make the link between the structural evolution of 
the melt and the variation of the order parameter with the pres- 
sure (^(P)). This will be done in a future publication. 
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